This semester I have been working in the Baldwin lab studying copepods (Aglaodiaptomus leptopus) from a bog in the northern Adirondacks near Sevey’s Point, NY. Past research shows that throughout the year, the copepods in this bog experience a change in pigment color. In the Spring they display red pigment and the color shifts to blue during the Summer months before returning to a red color in the Fall. My research explores what the driving force may be behind this pigment change and why. There are a few pre-cursors that are known to have effects on animal pigmentation which include diet, temperature changes, and UV light exposure. One that we are all familiar with is the effect that UV radiation has on our own skin. In the summer when UV radiation is most intense, it is common for our skin to become more tan due to increased melanin production. This fact is well-researched in humans, but with my research I wanted to examine if UV was also the driver of pigment change in these copepods.
Fig 1. Copepods collected on 10/20/22 at Seveys Bog
After collecting my copepod specimens from the bog, I anesthetized and photographed them under a microscope. I then imported those photos into Photoshop and recorded an average of the total body color of each specimen. In Photoshop, a color contains three values: L, a, and b which represent luminosity, the red vs. green scale, and the blue vs. yellow scale, respectively. After recording these initial values, I leave the zoopl. under different light treatments for 7 days and record how their pigment changes over the course of the week.
## 'data.frame': 706 obs. of 10 variables:
## $ image_name : chr "29Sept_D1a" "29Sept_D1b" "29Sept_D1c" "29Sept_D1d" ...
## $ L : int 64 63 62 66 64 77 74 66 68 80 ...
## $ a : int 2 2 2 6 2 3 6 3 2 2 ...
## $ b : int -1 -1 -1 8 -3 4 7 2 1 4 ...
## $ treatment : chr "D" "D" "D" "D" ...
## $ rep : int 1 1 1 1 1 1 1 1 1 1 ...
## $ before_after : chr "after" "after" "after" "after" ...
## $ date_collected : chr "9/29/22" "9/29/22" "9/29/22" "9/29/22" ...
## $ trial_duration_days: int 7 7 7 7 7 7 7 7 7 7 ...
## $ date_photo : chr "10/12/22" "10/12/22" "10/12/22" "10/12/22" ...
Workflow:
Plot > Guess relationships > Create model > Guess model assumptions > Run/interpret model > Replot to show key relationships
When I first brought the animals into the lab, I noticed anecdotally that they appeared to fade within a few days of being removed from their natural bog habitat. This led me to hypothesize that if I were to place the animals under total darkness (D) and a light treatment (NL) that contained little to no UV light, then I would be able to compare their luminosity values before and after a 7-day trial period. I used a fluorescent bulb as my “natural light” treatment to experiment with wavelengths in the visible light spectrum.
To examine the spread of data, I start by creating a data frame that includes the desired variables and converting character variables into factors.
luminosity <- DF[, c("L", "treatment", "before_after", "date_collected")] %>%
filter(DF$date_collected == "10/14/22")
luminosity$treatment <- as.factor(luminosity$treatment)
luminosity$before_after <- as.factor(luminosity$before_after)
luminosity$date_collected <- as.factor(luminosity$date_collected)
## L treatment before_after date_collected
## Min. :46.00 D :78 after :115 10/14/22:155
## 1st Qu.:57.00 NL:77 before: 40
## Median :61.00
## Mean :60.31
## 3rd Qu.:64.00
## Max. :75.00
Looking at the spread of luminosity (L) values, they are what I would expect. After working on color analysis all semester, I can confidently say that most L-values for my copepod images fall in the range of about 45 to 80.
A quick boxplot because displays the means and gives me a good idea about the variability within the treatments.
##### Guess relationships
For the darkness treatment (D), average luminosty before is about 54 and after is about 63.
For the natural light treatment (NL), average luminosty before is about 55 and after is about 62.
Mean luminosity before and after seem to be significantly different
Luminosity (L) values for darkness (D) and natural light (NL) treatments look to be fairly similar both before and after the treatments were applied, so I hypothesize that the treatments did not have a significantly different effect on luminosity.
##### Create model and guess assumptions
To check if this is statistically true, I will make a statistical model. I am choosing to use a generalized linear model because L-values are integer valued and are bound on a scale from 0-100. Therefore, a general linear model would not be appropriate because normal distribution cannot be assumed.
model_L <- glm(L ~ treatment * before_after, data = luminosity, family = poisson)
A few important things to note about my glm…
I have made L a function of treatment * before_after to test the interaction between these two variables. Later on when I run the statistics, I will be able to see if the treatment affects the luminosity differently when comparing before and after values.
I have decided to set the family equal to poisson because the L-scale on the y-axis is not continuous. As I mentioned before, the data are bounded by a scale, so the Poisson distribution accounts for this.
To ensure that the assumptions I am making about the distribution of my data are appropriate, I will use the autoplot function:
autoplot(model_L, smooth.color = NA)
The residuals appear to be evenly distributed and the data points follow a normal distribution according to the Q-Q plot, so I can confirm that my assumptions were okay and I can move on to running my model.
##### Run and interpret the statistical model
anova(model_L, test = "Chisq")
## Analysis of Deviance Table
##
## Model: poisson, link: log
##
## Response: L
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 154 72.726
## treatment 1 0.026 153 72.700 0.8711
## before_after 1 36.154 152 36.546 1.824e-09 ***
## treatment:before_after 1 0.785 151 35.761 0.3755
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
There are a few conclusions I can make based on this model:
The luminosity values between the D and NL treatments is not significantly different (p = 0.8711).
For both treatments, luminosity of the copepods increased significantly after 7 days (p < 0.05).
The light treatments did not affect the luminosity differently during the 7-day trial (p = 0.3755)
To visualize these findings and highlight the key relationships, I will create an interaction plot. The first step in this process is to create a separate data frame summarizing the mean values of luminosity for both treatments both before and after the 7-day period. Once I have those mean values isolated, I can use them to create the interaction plot.
Typically on an interaction plot, intercepting lines indicate that there is a significant interaction between two factors. For example, looking at this plot, it may appear that the darkness treatment causes the copepods to become more luminous over time than the natural light treatment, but as we confirmed with our statistical model, that is not the case. The overlapping of the error bars for the two treatments also indicates to us that the slopes of these lines could be variable and therefore we cannot say that one treatment significantly affects luminosity differently than the other.
luminosity does not change significantly.On 11/3/22, I collected copepod samples at the bog and photographed them 24 hours later to allow them to acclimate to the temperature in the lab. I analyzed these “before” images and will be referring to the L, a, and b values from this group as the control treatment. I then placed the copepods into three treatment groups: darkness (D), natural light (NL), and UV light (UV). After 7 days, I photographed these individuals and analyzed their color in order to compare them to the control.
##### Preliminary steps
To start, I will make a data frame with the variables I want to look at (L, a, b, and treatment).
Nov3 <- DF %>%
filter(date_collected == "11/3/22") %>%
select(L, a, b, treatment)
The “before” group is listed as “NA” for treatment, so I want to change that so it is listed as a “control” treatment. Then I want to convert the treatment variable into a factor
Nov3$treatment <- Nov3$treatment %>%
replace_na("control")
Nov3$treatment <- as.factor(Nov3$treatment)
##### Plot
The next step is to make a quick plot to look at the general spread of data. A boxplot will show the means of each treatment in relation to each other.
ggplot(Nov3, aes(x = treatment, y = L)) +
geom_boxplot() +
xlab("Treatment") +
ylab("Luminosity") +
theme_bw()
##### Guess relationships
Based on my work from earlier in the season, I would guess that luminosity increases significantly for the D and NL treatments, but does not for the UV treatment.
Similar to the previous test I ran, I used a generalized linear model with the family poisson and check my assumptions by looking at the residuals and Q-Q plot.
model_L <- glm(L ~ treatment, data = Nov3, family = poisson)
autoplot(model_L, smooth.color = NA)
All looks good, so I can move onto running my stats.
Since there is one categorical variable (treatment) and one integer variable (L), I will run a simple one-way ANOVA test to compare the three light treatment groups against the control group.
anova(model_L, test = "Chisq")
## Analysis of Deviance Table
##
## Model: poisson, link: log
##
## Response: L
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 238 62.275
## treatment 3 11.273 235 51.002 0.01034 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
The low p-value indicates that between treatments, there is significant difference in luminosity. To identify where this significant difference lies, I can perform a Tukey test.
modelL_tukey <- aov(L ~ treatment, data = Nov3)
TukeyHSD(modelL_tukey, conf.level=.95)
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = L ~ treatment, data = Nov3)
##
## $treatment
## diff lwr upr p adj
## D-control 5.500000 3.1338765 7.8661235 0.0000000
## NL-control 3.533333 1.1334867 5.9331800 0.0010151
## UV-control 2.161538 -0.1949861 4.5180630 0.0850771
## NL-D -1.966667 -3.5534886 -0.3798448 0.0082668
## UV-D -3.338462 -4.8589706 -1.8179525 0.0000002
## UV-NL -1.371795 -2.9442679 0.2006782 0.1112028
The copepods in the NL and D treatments became significantly more luminous. The UV copepods stayed darker and did not fade significantly from the original color.
##### Replot to show key relationships
ggplot(Nov3, aes(x = treatment, y = L, color = treatment)) +
geom_point(size = 3, alpha = 0.3) +
geom_point(data = mean_df_L, aes(x = treatment, y = mean_lum), shape = "diamond", size = 8, show.legend = FALSE) +
geom_hline(data = mean_control_L, aes(yintercept = mean), linetype = 2) +
theme_bw() +
xlab("Treatment") +
ylab("Luminosity (L)") +
coord_flip()
A scatterplot highlighting the mean L-values of each group effectively shows how the treatments differ compared to the mean.
a-value for color does not decrease significantly for this treatment.Another way of analyzing pigment is by looking at the a-value, which represents red-green scale where positive values represent red pigment and negative values represent green pigment. As I established earlier, the copepods fade significantly under darkness and visible light, but retain a darker shade under UV. In the Fall, the copepods at Sevey’s Bog are typically a red color, so I believe that if they are fading under D and NL treatments, they will be losing this red pigment and therefore the a-value will decrease. On the other hand, if UV triggers a retention of pigment to protect the copepods from harmful UV radiation, I would not expect the a-value to decrease.
##### Plot
ggplot(Nov3, aes(x = treatment, y = a)) +
geom_boxplot() +
xlab("Treatment") +
ylab("a-value") +
theme_bw()
Like the previous model I made for lumnosity, I used a generalized linear model, but this time with the the family gaussian because unlike the luminosity scale, the a-scale contains negative values which does not work with the Poisson distribution I used previously.
model_a <- glm(a ~ treatment, data = Nov3, family = gaussian)
autoplot(model_a, smooth.color = NA)
All looks good, so I can move onto running my stats.
##### Run and interpret the statistical model
anova(model_a, test = "F")
## Analysis of Deviance Table
##
## Model: gaussian, link: identity
##
## Response: a
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev F Pr(>F)
## NULL 238 1777.2
## treatment 3 396.7 235 1380.5 22.51 7.576e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
The low p-value for treatment indicates that there was a significant difference in the a-values compared to the control.
To find out exactly which treatments were different from one another, I am running another Tukey test.
model_a_tukey <- aov(a ~ treatment, data = Nov3)
TukeyHSD(model_a_tukey, conf.level=.95)
## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = a ~ treatment, data = Nov3)
##
## $treatment
## diff lwr upr p adj
## D-control -3.1133333 -4.691578 -1.5350883 0.0000041
## NL-control -3.5848485 -5.185587 -1.9841096 0.0000001
## UV-control -0.9974359 -2.569278 0.5744065 0.3571741
## NL-D -0.4715152 -1.529953 0.5869223 0.6572338
## UV-D 2.1158974 1.101692 3.1301031 0.0000010
## UV-NL 2.5874126 1.538546 3.6362791 0.0000000
From this, I can see more clearly that the a-values for the darkness and natural light treatments change significantly over the 7-day period (p-values < 0.05), whereas the a value for the UV treatment does not change significantly (p = 0.357).
##### Replot to show key relationships
Retention of dark pigments in copepods under UV implies that UV wavelengths stimulate photoprotection in these animals.
Although UV light did not cause the copepods to fully transition back to a blue color, future work could experiment with greater intensity and/or longer exposure of UV light treatment to see if this dramatic color shift can be induced artificially.
Generalized linear models are key when you have data that are not assumed to be normally distributed.
To import a second image, I had to change the “img-knitr” in the second image to something different.
I still need to figure out how to put horizontal lines into my scatterplot.
At first, I was confused about the interaction plot because the lines can still intersect even if the interaction is not significant.